Large amplitude behavior of the Grinfeld instability 
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Amplitude expansions are used to determine steady states of a semi-infinite solid subject to the 
Grinfeld instability in systems with a fixed (wave)length. We present two methods to obtain high- 
order weakly nonlinear results. Using the system size as a control parameter, we circumvent the 
problem that there is no instability threshold for an extended system in the absence of gravity. This 
way, the case without gravity becomes accessible to a weakly nonlinear treatment. The dependence 
of the branch structure of solution space on the level of gravity (or density difference) is exhibited. 
In the zero-gravity limit, we recover the solution branch obtained by Spencer and Meiron. A 
transition from a supercritical to a subcritical bifurcation is observed as gravity is increased or the 
nonhydrostatic stress is decreased at fixed gravity. At given values of the system parameters, we 
find a discrete, possibly infinite, set of solution branches. This is reminiscent of dendritic or eutectic 
growth, where similar solution sets exist, of which only a particular one is linearly stable. Despite 
the high order of our expansions, the approach is restricted to relatively small nondimensional 
amplitudes (< 0.2), a disadvantage we can overcome by a variational approach that will be discussed 
in a companion paper. At the critical point, we find that not only the first Landau coefficient is 
negative but all of them up to the highest amplitude order (15) we could compute so far. 
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I. INTRODUCTION 

In the processing of solids, it is often inevitable that 
they are submitted to elastic stresses. During solidifica- 
tion, for example, solids contract or expand according to 
their thermal treatment, which in a confined geometry 
will generally lead to the creation of strain [jl], 0, ||, ||, |j| 
and the formation of defects. Epitaxial growth, a tech- 
nique of utmost importance in the fabrication of novel 
electronic materials, necessarily produces strain in neigh- 
boring layers of different composition whenever there is a 
nonvanishing lattice mismatch of the substances involved 

10. 

The appearance of stresses is not restricted to environ- 
ments dealing with technical applications. A beautiful 
natural example is the Giant's causeway in Northern Ire- 
land, where cooling lava has turned into hexagonal basalt 
columns, most likely under the influence of stress. Very 
similar structures but on the mm scale instead of the m 
one, were obtained in the drying of starch j8|. It is there- 
fore interesting to study the influence of stress on the 
growth, melting, or evaporation of solids. 

As has been known for some time, the initially flat 
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interface of a solid subject to an elastic field tends to be- 
come unstable against morphological perturbations if an 
active transport mechanism is present. This phenomenon 
is often called the Grinfeld instability after its rediscov- 
ery by Grinfeld, who recognized its universal importance 
A description for a special system had been given 
earlier fiofl . 

The dynamics of the Grinfeld instability is energeti- 
cally driven. Therefore, the problem is variational and 
amplitude equations can be derived from an appropriate 
free-energy functional. Since the basic modes appearing 
as eigenfunctions in a linear stability analysis are plane 
waves for a planar interface pd| | , the natural approach is 
to decompose the interface in terms of a Fourier expan- 
sion and to determine the amplitudes of the modes by 
minimisation of the free energy. This is equivalent to the 
calculation of steady states from weakly nonlinear am- 
plitude equations with the amplitude expansion taken to 
one order less than the energy expansion. Consistency 
considerations as discussed in ]l2| apply here, too. 

As we shall see later, an approach that yields much 
faster convergence towards a steady-state solution can be 
based on an expansion into a different set of generalized 
functions in the spirit of a variational approach. This 
is the multicycloid method to be presented in a second 
paper p3[ , henceforth referred to as II. Its utility goes 
farther than improving convergence. Certain interfaces 
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with overhangs are tractable within that approach, which 
extends its applicability beyond that of the amplitude 
method. 

In the model, there is competition of at least two en- 
ergetic terms: the strain energy stored in the bulk of the 
stressed material which tends to destabilize the surface 
at any given wave number, and the surface energy which 
is stabilizing. In many systems, more physics is added by 
the presence of potential energy in a gravitational field 
which acts stabilizing as well, if the density stratification 
corresponds to high densities (solid) below the interface 
and low densities (melt) above. 

The main goal of this paper is to find steady state 
solutions to the dynamical equations arising from the 
Grinfeld instability. Transport can be, for example, via 
"melting-crystallization" waves, when there is contact 
with the melt, constituting a particle reservoir, or by 
surface diffusion, for the case of a solid in vacuum. In 
the first situation, the dynamics is nonconservative, in 
the second, it is conservative. However, kinetic aspects 
of the transport mechanism do not play a role in the final 
steady state. Therefore, our results will be applicable to 
both types of systems. 

What is known about steady states for the Grinfeld 
instability so far? Spencer and Meiron |Q have shown 
by time dependent simulations using a highly accurate 
spectral code that for a solid submitted to stress in the 
absence of gravity there is a branch of stable steady- 
state solutions only for a very limited set of wavelengths. 
Moreover, they showed that the generic behavior for large 
wavelengths is the evolution of cusp singularities. These 
singularities occur for interface shapes which are very 
close to the shapes of certain cycloids, a result that is 
suggested also by the analytic work of Chiu and Gao 

We will present results on the solution space of steady 
states including gravity. As we will discuss, there appears 
to be a ladder of solution branches reminiscent of the sit- 
uation in dendritic and eutectic growth. In the former 
case, there exists a discrete infinite set of needle solu- 
tions for any given undercooling, of which however only 
(the fastest) one is stable. In the latter whole set 

of solution branches exists in a certain wavelength range 
for given temperature gradient and pulling velocity. In 
the present approach, different solution branches corre- 
spond to different wavelengths of the pattern. Similar 
morphology changes occur between branches as in eutec- 
tic growth p6| . 

Not including a full stability analysis using Floquet- 
Bloch theory, our approach allows us to make stability 
statements only for periodic perturbations. With respect 
to these perturbations, there is, in the absence of gravity, 
one branch that is stable along its full length, which is the 
one found already by Spencer and Meiron. With grav- 
ity, there is a critical wave number, below which small 
amplitude solutions become unstable. Only a part of 
the branch at finite amplitude is stable rendering the bi- 
furcation from the planar solution subcritical in a sense 



specified more precisely below. 

In principle, knowing the coefficients of the mode ex- 
pansion, we could extend our approach to time dependent 
calculations. This will however not be pursued here. 

The paper is organized as follows: Section II describes 
the basic method and gives expressions for the relevant 
energies of the system. In Sec. HI, we rescale the prob- 
lem for nondimensionalization and then formulate it more 
technically, adapted to the Fourier expansion. The next 
section, Sec. IV, compares the results of two methods to 



obtain the minima of the energy functional. One of them 
solves all the equations but one analytically, with the 
second, we have to solve the final set of nonlinear equa- 
tions for the amplitudes by a numerical scheme (Newton- 
Raphson method). In Sec. |y|, we add gravity leading 
to the understanding of the small-amplitude behavior 
down to the critical stress case described by Nozieres [JH] . 
Gravity plays an essential role for the solid-liquid inter- 
face close enough to the instability threshold |17|] . 



Finally, in section VII, we will summarize the conclu- 



sions that can be drawn from our high-order expansions. 



II. THERMODYNAMIC AL CONSIDERATIONS 

In general, in order to find steady solutions, we have 
to equilibrate the contributions of elastic, surface and 
gravitational energy terms. Basically we have to solve 



5E = 6 \E, 



clastic 



E 



surface 



E, 



gravity 



] = o 



(1) 



where E stands for energy or free energy and the sub- 
scripts describe the origin of the respective energy. Next, 
we have to verify the stability of our solutions. 

One way to express variations of the surface is to 
rewrite its position in terms of some amplitudes on re- 
lated to a given system of functions and to vary this dis- 
crete set of amplitudes. In the first part of this article, 
we will expand the interface shape into a Fourier series, 
assuming periodicity. It can be restricted to even (i.e., 
cosine) modes [Eq. (|2])] because of the reflection symme- 
try of the system about a; = 0: 
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CO 

E 



a n cos (nkx) 



(2) 



In a subsequent paper, we will consider an expansion 
in terms of multicycloids, which will prove to constitute 
a more appropriate tool for capturing cusp singularities 
at leading order. 

Our model is a two-dimensional system possessing an 
interface located about z — so that the mean interface 
position is zero. This condition is automatically satisfied 
if z is given by equation (|^). The uniaxial prestress cto is 
taken parallel to the direction of the abscissa, x, and the 
elastic domain is reaching downwards to minus infinity 
(whence our freedom of choice of the average interface 
position). 
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A. Strain Energy 

For linearly elastic solids, the strain energy density w 
is given by 



1 



1 



w 



-SijkiOijOki , 



(3) 



where <7y is the stress, ity the strain tensor and Siju is 
the compliance moduli tensor. As we assume an isotropic 
system and plane strain (i.e., we restrict ourselves to two- 
dimensional deformations), Syjy becomes simply 



S, 



1 



ijkl 



Ey 



(SikSji 



v8ij8ki) 



(4) 



energy of an infinitesimal slice, considering that the cen- 
ter of gravity is at £(x)/2, and by integrating then. Wc 
obtain 



9&p 



gravity 



(7) 



where Ap stands for the density contrast and g is the 
gravitational acceleration. 



III. NONDIMENSIONALIZATION AND THE 
DERIVATION OF HIGHER-ORDER EQUATIONS 

A. Nondimensionalization 



(see for instance [Q). Herein, v is Poisson's ratio, Ey is 
Young's modulus, and summation over repeated indices 
is implied. The total amount of elastic energy stored is 
obtained by 



E, 



clastic 



2Ey J-p 

(l + v)v 
E Y 



(JijOij dr 



det a dr , 



(5) 



where the domain T> is bounded by the surface and ex- 
tends to minus infinity in the z direction, whereas it can 
be restricted to one period of the basic mode in the x 
direction. It turns out that the determinant term is zero 
at each order of the mode expansion. Therefore, from 
here on, it will not be explicitly written anymore. The 
detailed method of calculating the strain energy will be 
given in section [II. 



B. Surface Tension 

Surface tension is assumed to be isotropic. The cor- 
responding total energy of the surface then is the length 
of the surface multiplied by the surface tension 7. If we 
take the zero of this energy to correspond to the planar 
interface, -E su rface of a deflected interface is: 



-^surface 7 



(Vl + C'(20 2 - l) dx 



(6) 



Here, as in Eq. (^) , we have to take a period of the surface 
for the integration interval B. 



C. Potential Energy 

Finally, the potential energy density, if referred to the 
planar interface as base state, becomes proportional to 
the mean square amplitude of the surface shape. This 
can be verified very easily by writing down the potential 



As long as gravity is not considered, all physical pa- 
rameters can be scaled out of the equations of motion 
for the Grinfeld instability, which renders the problem 
scale invariant. The length scale li to be used in this 
pro cedure, is, up to a factor, equal to the Griffith length 
|I9| , a measure for the competition between capillarity 
and elasticity. Setting wq = Cq(1 — v 2 )/2Ey-, which is 
the elastic energy density of the prestressed planar state, 
it reads lEl 



2w 



(8) 



When gravity is present, a second length scale I2 ap- 
pears, given by 



h = 



3Ap 



(9) 



An important parameter of the nondimensionalized equa- 
tions is the ratio of these length scales 



h2 := ^ 



h 

2l 2 



(10) 



Note that in the absence of gravity, this ratio becomes 
zero as the gravity length I2 diverges. 

Using these quantities, we nondimensionalize equa- 
tions (|), (§) and (0), carrying out a formal transfor- 
mation x — > l\x, z — > l\z, k — > Zj~ 1 k, £ — > Zi£ and 
cry — > <7o<7ij. Thus, in equation (||), we get a Cg and 
an additional h coming from the integration in the z di- 
rection - resulting in a new prefactor 7/2. Equation 
just keeps the prefactor 7, and Eq. (Q) can now be written 
with the prefactor 7Z12/2, accounting for an l\ 2 coming 
from the integrand. The x integration produces a pre- 
factor common to all integrals that can be dropped. 

Dropping the common prefactor 7 as well, our nondi- 
mensional equation reads 



5x 



dx 



= . 



(11) 
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This equation is still exact, of course only within the 
framework of our basic assumptions (no cross-effects be- 
tween stresses, capillarity and gravity, i.e., no capillary 
overpressure, stress tensor divergence free, etc., see [po|). 
The following calculations are however done assuming 
that the square root can be expanded in powers of (^'(x) 2 . 
Essentially, this restricts the validity of results to inter- 
faces with |£'(a;)| < 1- If we do not want to consider 
gravity, we can just set ?i2 = 0. It should be kept in 
mind that from now on we always consider nondimen- 
sional quantities. 



Eq. (g). We introduce an expansion parameter e and 
assume a n cx e n . As is well-known, in two dimensions 
calculations can be facilitated by introduction of an Airy 
stress function x( x i z )> connected with the stresses via 



Fx 

dz 2 ' 



d 2 X 
dx 2 1 



d 2 X 
dxdz ' 



(15) 



i.e., our stress function refers only to the deviations from 
the constant prestress. 

Using the boundary conditions at z = — oo and solving 
for the z dependence, we obtain 



B. Expansion of the stress components into Fourier 
modes 



To compute the elastic energy, we obviously first have 
to solve the elastic problem, i.e., the Lame equations, 
with appropriate boundary conditions. Next the stress 
tensor is expanded into a trigonometric series, exploiting 
periodicity in the x direction. The elastic partial differ- 
ential equations can then be reduced to ordinary differ- 
ential equations for the z dependent coefficients. After 
solving for the z dependence analytically, we can express 
the stress tensor in terms of expansion coefficients for 
arbitrary surface shape C,{x). This is a useful step, ren- 
dcring feasible the minimization procedure expressed by 
Eq. Q. 

Starting with an undisturbed plane half space, the 
nondimcnsional prestress takes the form: 



So = 



1 





(12) 



3G 



and we impose the boundary condition that for z - 
the total stress tensor approaches this quantity. 

The changes of the stress tensor due to a surface cor- 
rugation will result in a relaxational stress £ a . As the 
surface of the solid is free - the solid is in contact with a 
liquid - the total stress tensor £ = So + S has to fulfill 
the boundary condition 



(So 



Sn = —pen = 



(13) 



at the interface, where we have, for convenience, set the 
pressure pe of the liquid equal to zero. This is possi- 
ble, because the dynamics of the system is unaffected 
by an additive constant pressure | pp| . It is here, where 
we neglect capillary overpressure, which would impose a 
(small) jump discontinuity of the component S nra . 

Calling the components of S Q a xx , a xz = a zx and a zz 
and leaving out a common prefactor in the normal vector, 
Eq. (pi) reads 



(1 + <Jxx) C( x ) + °xz = 
-&xz C(x) + G zz = . 



(14) 



where stress components have to be taken at the inter- 
face and £(x) is the free surface, expressed according to 



n=o 



(a n + b n z) cos (nkx)e~ 



(16) 



where we anticipated the leading e-dependence of the co- 
efficients. Now, a n and b n have to be computed by use 
of the boundary conditions (|l4|). To this end, we first 
expand a n and b n : 



oo 

E ^ 

m=0 



n h — \ " c 2m h 



2)n. 



(17) 



m=0 



Note that odd expansion coefficients disappear, which 
is proven in appendix (|^). Furthermore, equations Jl4| ) 
are given at the interface. So, the exponentials in the 
stress components have to be expanded as power series in 
z, and then z has to be replaced by the interface position 

ax). 

After this has been done, equations (Jl4|) will be two 
power series in e, and for every order e™, there exists a 
set of equations provided by the contributions of different 
cosine and sine modes, respectively, which have to be set 
equal to zero separately p3[. 



C. Second and fourth-order results 

To second order in the amplitude, the sum of the three 
energies is 



AE= (-2k + k 2 + l 12 ) -i 
v ' 4 



(18) 



We write AE from now on as a reminder that we have 
taken the planar interface as the reference state with en- 
ergy zero. 

From Eq. ([l8|), we can recover two well-known results: 

• When neglecting gravity we get marginal wave 
numbers k — {0, 2}. 

• Under gravity the critical wave number is k = 1, 
and the critical length ratio Z12 = 1. 

More generally, we can conclude that for arbitrary Z12 
in the interval [0, 1] the band of unstable wave numbers 
ranges from 1 — VI — I12 to l + — Z 12 - If Z12 > 1, which 
corresponds to low stress or high gravity, the flat surface 
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is linearly stable. If we now add fourth-order terms we 
obtain 

where we have set A n = fca„ for short. The term hav- 
ing the prefactor 1/k is the elastic contribution to the 
free-energy change, the term with li 2 the gravity one, 
while the remaining three terms describe the surface en- 
ergy contribution. Minimizing this result with respect 
to A 2l the weakly nonlinear expansion of Nozieres Jll[] 
is recovered. The main result is that the instability is 
subcritical: Even for li 2 slightly larger than one a finite 
perturbation may trigger the instability. 

IV. EXPANSION METHOD 
A. Principles 

Having calculated the change of energy to a desired 
order 2n, we can now devise methods for solving the sys- 
tem 

A,:=^=0; i = l...n. (20) 

A scheme that goes a long way analytically is what we 
call the expansion method. Here we expand, besides the 



coefficients a n , b n , also the amplitudes a.% . . . a n , either 
in powers of e or, more straight-forwardly, in powers of 

00 

a m = a™ ^ c m ,j oP\ ■ ( 21 ) 

j=o 

If we calculate AE using the first n amplitudes, then 
it seems at first sight that we should actually expand 
Aj up to a 2 ™ -1 . But then we would have to take into 
account terms up to order a 1 n ~^ 3 ~ in the expression for 
AE before taking the derivative (which reduces the order 

by j). 

However, the structure of the equations is such that 
obtaining a n to order e n from setting A„ = 0, we get 
ttn—i to order e n+1 from the equation A„_i = (in which 
each a n is multiplied by at least a factor ai), then a n _2 
to order e n+2 , and so on, until the first equation (Ai = 0) 
yields ot\ accurate to order e 2n_1 , as desired. 

The coefficients of the powers of o.\ in A2 . . . A„ pro- 
vide a linear system of equations for the calculation of all 
c m ,j with m + j < n. Interestingly, similar to the Airy 
coefficients, all c mj - where j is odd are zero, a formal 
fact allowing us to considerably increase the speed of a 
calculation. 

Finally, Ai becomes a function of ol\ and k only. So we 
can numerically solve Ai = for any given wave number. 

In order to provide a flavor of the shape of these terms, 
we give the general solution for n = 4, using the short- 
hand expressions from equation ( |l9|) again: 



a\ 1 -294fc 3 + 234fc 4 - 94kl 12 + 121k 2 l 12 + 10f? 2 - 252fc 2 4 

/i 2 -4fc + 4fc 2 3 {l 12 - Ak + Ak 2 ) 2 (l 12 -6k + 9k 2 ) ^ 

1 -60fc 2 - 80k - 12Zi 2 + 12k 3 + ikl 12 3 
8 (;i2-4fc + 4fc 2 )(/i2-6fc + 9fc 2 ) ai 

1 -288/c 6 + 644fc 3 /i2 + AAkl\ 2 + 2112fc 5 - 72k 2 l 12 - 672fe 3 + 13fc 2 Z 2 2 + 2Z 3 2 - 960fc 4 - 52fc 4 / 

2 (l 12 -4k + 4fc 2 ) 2 (l 12 - 6fc + 9k 2 ) (Z12 - 8fc + 16A; 2 ) 

I 



We recognize a restriction of our method: because of 
the denominators originating from the solution of our 
system for a m j the functions representing a n will contain 
a number of poles, increasing with order. Their locations 
depend on the value of Z12. Generally, the positions of 
poles are solutions of 

h 2 -2jk + ] 2 k 2 = 0; j = 2...n, (23) 

that is 

Kpole = : ■ (24) 



For the no-gravity case li 2 = 0, we find poles at 2 /j 
(j = 2 . . . n) . So the method will be reliable for all k 
values sufficiently far above k = (1 + VI ~ ^i2)/2 (unless 
the solution of the final nonlinear equation leads to a new 
singularity). If li 2 > the zeros move to the left, thereby 
increasing the range of validity of our method. 

Hence, at the critical wave number k = 1 of the model 
with gravity, the method is usable because the largest 
k value with a pole is located at k = V 2 f° r '12 = 1- 
In Fig. 0, we have plotted the positions of the marginal 
wave numbers and the first poles, depending on l± 2 . The 
range of k values, in which the expansion method should 
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FIG. 1: Location of the marginal wave numbers (bold) and 
the first pole coming from the denominators of the amplitude 
solutions. 



work properly, is bounded from below by the upper part 
of the thin line. 

Another point we have to stress concerning the valid- 
ity of solutions results from the necessity of expanding 
the square root in the arc length formula [Eq. @], giv- 
ing the surface energy. This expansion does not converge 
for C( x ) 2 > 1- As long as the integral of the resulting 
divergent series is still convergent, the total energy will 
not be strongly influenced by a nonconvergence of the in- 
tegrand. Therefore, a violation of this criterion within a 
small x interval may not immediately render the calcula- 
tion useless. We find that in practice numerical solutions 
start to diverge when the criterion 



(25) 



is not satisfied. The angular brackets denote a spatial 
average over one periodicity unit. 

A third factor which might become problematic will 
be pointed out briefly: The solution of Ai = 0, which 
has to be carried out numerically has to deal with very 
large natural numbers appearing in the problem when 
we go to high order. For a reliable solution the computer 
algebra system (CAS) Maple is an appropriate tool for 
two reasons: The engine works with fractions of natural 
numbers as long as possible, thereby avoiding truncation 
errors, and for the genuine floating point operations we 
can adjust the accuracy as necessary, which is essential 
for the results to follow. 



B. Solutions with gravity neglected 



The method introduced in section IV A will now be 



applied to the simplest variant of the model where only 
elasticity and surface tension are present. Therefore we 
set l\2 = 0. A first check is a comparison with the nu- 
merical result of Spencer and Meiron |L4| |. 

Figure ^| shows how the method approaches the numer- 
ical results. The shapes corresponding to each data point 



FIG. 2: Approximation of the numerical solution by the ex- 
pansion method with accuracy n = 2 to n = 8. The bold 
line is the (digitized and scaled) numerical solution from fL4| . 
Numbers in circles correspond to curves in Fig. ^. 




FIG. 3: Shapes of some particular solutions of the no-gravit 
model calculated using eight amplitudes, according to Fig 
Curves are scaled in the x direction to map them onto an 
interval of length 2-k. 



grow from a shallow cosine to a sharply grooved struc- 
ture, converging to the cycloid-like behavior proposed by 
Chiu and Gao M, see Fig. ||. 

With increasing order the solutions get closer to the 
numerical branch. Remember that every solution point 
corresponds to a whole shape. The shapes corresponding 
to marked points in Fig. ^ are shown in Fig. [s| The more 
the solution approaches a cusped structure, the more in- 
accurate it becomes due to the lack of higher modes nec- 
essary to describe the cusps correctly. A sign of decreas- 
ing accuracy is the existence of nonphysical undulations 
in the shapes, as observable in the outer parts of curve 



Moreover, we point out that all solutions are stable 
which can be checked by calculating the second deriva- 
tives with respect to the amplitudes contained. 
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C. Secondary solutions in the absence of gravity 

Besides the solution which has already been found nu- 
merically [Q, our method finds several types of sec- 
ondary solutions whose physical relevance has to be dis- 
cussed. 



1. Solutions at large k values 

First, there appears to be an infinite number of sub- 
critical "bubbles" at increasing wave numbers. To con- 
vey an idea of what they look like we present the first 
of these branches at different orders of the calculation in 
Fig. ||. Obviously, to the order of accuracy available to 
us, namely 8, it is impossible to decide with certainty 
what will be the behavior of these solutions at infinite 
order. 

At order two there is not a complete bubble but a 
branch decaying to zero at infinity. This can be seen by 
inspecting the second-order energy, Eq. (19), replacing 
Ai by the minimized A2 calculated in Eq. (p8|), taking 
the derivative with respect to A\ and setting l\2 = 0. 
Apart from the trivial solution A\ = we find 



A\ 



2-3k + k 2 
3fc 2 - 27/c + 40' 



which goes to 8/3 at large k. Calculating the mean square 
amplitude a is connected with a division by k, hence the 
decay to zero. 

At order three, this branch bends up again, diverging 
at k = 22.45 with no real solutions for higher k and at 
order four it moves to lower k (this case is shown in part 
A of Fig. |^) , and an additional branch (not shown) at 
higher k emerges, also bending up at order five and so 
on. 

Let us have a look at the energy formulas. As we have 
expanded them in powers of A±, they are polynomes. 
Consequently, the large amplitude behavior is dominated 
by the coefficient of the largest power of A\. Moreover, it 
turns out that for k > 2, the lower powers in the solutions 
do not produce additional extrema. 

So it is easy to decide about the stability of the solu- 
tions. As we know the flat surface to be linearly stable 
for k > 2, the next solution towards higher amplitudes 
at any fixed k should be unstable and the next one, if 
it exists, stable again, because energy minima and max- 
ima must alternate as a function of the amplitude. So all 
branches in Fig. |] are stable above the symbols, since 
the energy is minimal there. 

An interesting conclusion is suggested by figure |]D. 
The onset point of the first additional branch seems to 
converge to k = 2. As with no method available another 
stable solution at this point is found, the branch might 
become infinitely steep at high order, which would mean 
that the exact energy landscape becomes independent of 
amplitude at k = 2. 
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FIG. 4: Additional solutions at different orders. In each case, 
not all additional branches are shown. The symbols denote 
the change of stability. For comparison, we have included the 
numerical solution from jL4| in each plot. A: First additional 
branch at order four. The second, decaying one starts above 
k « 60 here. B: First two additional branches at order six. 
The third, decaying one starts above k ~ 329 in this approx- 
imation. C: Only the first additional branch at order eight. 
There are three more solutions not shown in this diagram. 
D: Onset points of the first additional stable branch to orders 
four through eight. 




FIG. 5: An example shape of the additional solution at the 
onset of the stable branch in Fig. [l|c, stretched to an interval 
2tt. 



A less spectacular interpretation would be that the ad- 
ditional branches just get very closely spaced without 
collapsing onto single lines. 

Figure |s| gives an interface shape for the first additional 
branch, showing explicitly that interface morphologies 
from this branch are significantly different from those of 
the original branch (compare with Fig. ||) . While this so- 
lution is stable with respect to amplitude perturbations 
of the same periodicity, we expect the whole branch to 
be unstable with respect to large wavelength perturba- 
tions. A Floquet-Fsloch analysis would presumably show 
that each new branch has one more unstable eigenmode, 
a situation that is also encountered in dendritic growth 



8 



when analysing steady-state needle crystal solutions. 



2. Small k values 

First of all, we should notice that for a curve £(x) given 
by its Fourier series, which is a solution, there exists an 
infinite number of alternative representations 



C(x) = ctj cos ( nj 



n € N. 



(26) 



So, if we for instance choose n = 2, then we should expect 
a "ghost branch" starting at k = 1 instead of 2 which 
contains only modes 2, 4, 6 . . . . As our expansion method 
requires the first amplitude, A\, to be nonzero, we did 
not expect those solutions to be found by it. 

Yet, notwithstanding these considerations, our algo- 
rithm finds solution branches at k = 1, k = 2 /3, etc., 
with a shape of the mean square amplitude plot similar 
but not identical to that of the branch at k = 2. To 
understand this, let us consider the case k ~ 1. Numeri- 
cally, Oil can become very small while still being nonzero. 
Then u-i will be finite - the first denominator in Eq. ( p2| ) 
is close to zero -, will be small again, because the de- 
nominator is 0{a\ ) and the numerator 0{a\ ), will be 
finite, and so on. Hence, we just obtain the ghost branch 
at n — 2 from Eq. ( prf ) . Clearly, the numerical quality of 
approximation will not be good in these circumstances. 
At larger amplitudes, the ghost branches are not exact 
images of their original anymore, because the Fourier ex- 
pansion will contain increasingly large contributions from 
the "forbidden" modes cos (Ikx/n), I ^ nj. 

By use of the numerical scheme described in section [vj, 
one avoids the zeros of the denominators, stemming from 
the expansion of the amplitudes in powers of a±. Instead, 
the denominators will be renormalized by the amplitudes 
themselves. Hence, it is possible to reproduce the exact 
image of the k — 2 branch as a ghost at k = 1 and at 
lower k values. Moreover, one may decide whether there 
are new truly different solution branches near the ghosts. 
It turns out that this is not the case, i.e., the only true 
solution branch in the lower vicinity of k = 1 is the n = 2 
ghost according to Eq. (p6[). 

We conclude that in the case without gravity, the 
Spencer-Meiron branch is presumably the only one that 
is stable in finite systems. All additional branches are 
either "ghost images" of this one or are unstable. 

It must be kept in mind, however, that the existence of 
steady-state solutions, even though they may be unsta- 
ble, could play a decisive role in dynamics. An example is 
given by the Kuramoto-Sivashinski equation [2lf| , where 
there exists an infinite set of unstable steady-state solu- 
tions for which it has been attempted to connect them 
with spatio-temporal chaos p2]| . 




0.40 
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0.80 
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1.00 



FIG. 6: n — 2, 3 and 4 ghost branches, emerging at k = 1, 
2 /3 and !/2, respectively. 



V. THE CASE WITH GRAVITY 

A. The critical case 

The classical linear analysis shows that at a critical 
stress 



E 



■V A P91 



a single wave number 



k c . — 



lApg 
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becomes neutral. In our reduced units, this critical point 
is at 



1 



1 . 



Nozieres calculated the change of the free energy up to 
order af |ll[] at this point. The negative sign of the 
coefficient of the a\ term lead to the conclusion that the 
Grinfeld instability is subcritical. As a first example, 
we calculate this quantity up to a\ 6 using our method. 
As usual, the prefactor 7 which is also contained in the 
elastic energy due to the rescaling has been divided off, 
and so, at k = 1, the change of the free energy is 



AE 



43 4 
"64 ai 



11661 
1024 



3446989 



O'l 



317790953065 



12288 
154676894913120211 



■ ai 



37748736 543581798400 
81526198232945992800847 14 

7827577896960000 ai 
2229374504089574965742093912371 



16 



5523138964094976000000 ' Ul ' ^ 

The prefactors are negative at all orders, and their abso- 
lute value increases exponentially. This behavior shows 
that a restabilization does not occur, at least not at the 
basic wave number k — 1. For the elastic energy can be 
reduced indefinitely by increasing the amplitude of the 
system p3. 
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B. Analysis in the full £12 range 

We will now consider nonzero values of l±2- Solution 
branches will start at wave numbers smaller than 2 as 
has already been pointed out in discussing Fig. |l]. 

In Fig. we show the result of the highest-order expan- 
sion available so far, including eight amplitudes. In doing 
so, we restrict ourselves to the branches emerging at up- 
per marginal wave numbers. There are more solutions 
for < k < 1. 




FIG. 7: Solution branches for different I12 values (small num- 
bers), obtained using the expansion method. Solutions above 
the dashed line are probably less reliable because at least in 
the no-gravity case the analytical result starts to consider- 
ably differ from the numerics (bold again). Other symbols 
are explained in the text. 



We will now consider the weakly nonlinear analysis for 
arbitrary wave numbers taking gravity into account. The 
small-amplitude behavior can be entirely derived from 
the simplest weakly nonlinear model, namely Eq. (19). 
For steady solutions we must take AE from this equation 
and require 



dAE 



0, 



dAE 



0. 



which leads to the condition 



2A\k 



4fc(l - k) -I 



12 



(28) 



Inserting this into the second equation results in a rather 
lengthy expression from which at any (fc, Z12) the A\ part 
of steady solutions may be obtained. A2 is then calcu- 
lated by equation (28), and the quantity plotted in Fig. 
is the mean square amplitude of the resulting pattern, 
which is defined by 



In the special case considered, we simply have 

W S = Z=M±5. (30) 



Now we would like to discuss the stability at the onset of 
the branches. Remember that the flat interface is stable 
against wave numbers larger than the upper marginal 
one. The nontrivial steady branch emerges right at the 
marginal wave number, and for sufficiently low values 
of 1 12 it continues at finite amplitudes towards smaller k 
values. For I12 — and small amplitudes we find a square 
root behavior. 

This branch pattern, which we may call supercritical 
(although the wave number k is not what one would usu- 
ally talk about as an order parameter) changes to sub- 
critical at a value which may be determined by consider- 
ing the branch emerging uncurved from the bifurcation 
point. Hence, to obtain this tricritical point, we just have 
to calculate the parameters, where the second derivative 
of AE with respect to the amplitude A\ vanishes at the 
onset of the nontrivial branch, i.e., we have to solve the 
system: 



1 + Vl-'l2, 



d 2 AE 
8A\ 







The solution is 



lll.t — — 



148 20 



9 



9 



57 « 0.333 



13 1 



kt = — - 1.817. 



(31a) 



(31b) 



Consequently, all branches emerging left of k t (bullet 
symbol) in Fig. ^ are unstable at the bifurcation point. 

Another remark has to be made: This simple weakly 
nonlinear model is accurate for amplitudes up to w 0.2. 
So, solutions obtained at larger amplitudes should not 
be trusted. In II, we will present a method for getting 
reliable solutions at larger amplitudes. 

Now we consider the stability of branches emerging 
to the left of k t . It can be examined by differentiating 
AE(Ai, k) by A\ twice and computing its sign at all the 
data points. 

For a parameter range up to l\% = 0.6 this yields onset 
points of stability, which are located, where the inclina- 
tion of our solution branches changes sign (open circles 
in Fig. 0). At even higher l\ 2 values, a possible onset of 
stability cannot be calculated anymore as it would lie in 
amplitude regions which are not correctly treated by our 
expansion. 

As a first conclusion we would like to point out that 
the results in Fig. (?] lead to the analog of the Spencer- 
Meiron branch in a system with gravity, showing it to be 
stable for small gravity (with respect to fixed-periodicity 
perturbations) at small amplitudes, whereas for larger 
gravity stability begins only at a finite amplitude. 

Figure || shows the corresponding branches emerging 
from the lower marginal wave numbers, for different l\ 2 
again. At the onset they are all unstable, and the large- 
amplitude behavior is difficult to assess because there is 
a rapid change in the amplitudes at almost the same k 
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value. The shape of solution branches is very sensitive to 
changes of order, and for the larger Z 12 parameters our 
method stops finding solutions at all. 
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FIG. 8: Branches emerging from lower marginal wave num- 
bers at different £12, to highest accuracy available. Comments 
see text. 




which are not shown here. This behavior can be ex- 
plained by regarding Fig. [|C again: at k = 1.9 the 
unstable solution branch diverges, and the method pre- 
sented in this section obviously produces a nonphysical 
connection between the two branches. 



VI. NUMERICAL INVESTIGATION 



VII. SUMMARY 



Having calculated the change of energy due to a cor- 
rugation up to some amplitude we can apply an alterna- 
tive method to calculate our sets of amplitudes. This is 
to take the energy formula and calculate its derivatives 
[Eq. (p0|)] without any expansion. Then hx k and try to 
find a solution vector A = (A±, A2, . . . ) numerically. 

When applying this method, a number of nonphysi- 
cal solution branches emerge. This happens because of 
the restricted set of equations which is used: if we take, 
for instance, the energy to sixth order, then it contains 
the first three amplitudes, and there will be a term like 
A\A2, but not A\A2- So, artificial intersection points 
in the parameter space {A\, A2, A3} are produced which 
are smeared out and replaced by other artificial solutions 
at higher orders. 

We therefore have to start searching at k slightly below 
2, thereby restricting the solution space to small absolute 
values. This ensures that the solution found is on the 
branch already known from (Ml and section IV B. By 



slowly varying k and taking the previous solution vector 
as a guess for the new one we can scan a whole k interval, 
until the method breaks down. The procedure can be 
programmed with a few lines of Maple code, although 
the values of the different solution amplitudes differ by 
several orders of magnitude and the system of equations 
is highly nonlinear. 

We have calculated solutions again for the simplest 
model, at I12 = 0. The results are shown in Fig. ||. 

The solutions possess a good qualitative agreement 
with the numerics up to k slightly below 1.9 where, at 
least for order 3 and higher, the numerics suddenly ter- 
minates by jumping to one of the nonphysical branches 



We have extended the weakly nonlinear analysis of the 
Grinfeld instability to higher order in the amplitude than 
previous approaches p| . 

An interesting fact we found here is that at the crit- 
ical point (k c , l\2, c — 1) all studied orders in the am- 
plitude expansion (up to the 15th) showed subcriticality 
[see Eq. (p7|), which gives the energy at threshold to 16th 
order]. We conjecture that this phenomenon continues 
to infinite order. This feature is rarely encountered in 
nonlinear systems and interesting in itself. 

For the Grinfeld instability in the absence of gravity, 
this is the first time such an expansion is performed. 
Since in that case, there is no instability threshold for 
the infinite system, the system size/wavelength or else 
the wave number may be taken as a control parameter. 
We recover the Spencer-Meiron result - a supercritically 
bifurcating branch of steady-state solutions, which end 
in a cusp singularity (that we cannot quite reach with 
our approach). 

As gravity (i.e., I12) is increased, the bifurcation moves 
to smaller wave numbers and changes from supercritical 
to subcritical. Stable finite-amplitude solutions are found 
for intermediate values of I12 but not for all of them. For 
example, between I12 = 0.6 and l\2 = 1.0, there are no 
stable solutions with amplitudes smaller than 0.2, and 
for larger amplitudes, the methods presented here are not 
really useful. In II, we will be able to make more precise 
statements about the behavior at larger amplitudes, as 
the variational approach used there is more powerful. 

There is a second bifurcation in wave number space at 
small wave numbers, where the dispersion relation crosses 
zero again (for nonvanishing gravity, the k = mode 
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must always be linearly stable). All branches emerging 
there are however unstable. Also we cannot interpret this 
bifurcation as being controlled by the system size, be- 
cause small wave numbers correspond to large systems, 
which already contain unstable modes of smaller wave- 
length. So this second branch is of no physical relevance. 

With the mode expansion, one can find more branches 
of solutions, which do not bifurcate from the planar state 
but lead a detached existence. They might have relevance 
in very small systems which are linearly stable but forced 
into a new state by a strong perturbation. However, we 
cannot make clear-cut statements about the stability of 
these states, since the expansion converges very slowly. 
As we have pointed out, even if they are unstable, these 
states may by their very existence play an important role 
in the dynamics. 

To summarize, this approach is useful for identifying 
possible stable steady states (which are only metastablc 
in an extended system), but despite the high order ac- 
cessible in the expansion using computer algebra, it is 
limited in the size of the amplitudes about which reli- 
able assertions can be made. The variational approach 
to be discussed in II goes beyond this limitation and is 
even capable of producing cusped solutions analytically. 
Moreover, it is not restricted to interfaces without over- 
hangs. 



variance. The derivative is then given by: 
2 

e n nA n sin (n£) 

71=1 



C'(0 



(A2) 



The simplest term is gravity, and the nondimensionalized 
and transformed integral (|7j) is given by 



(A3) 



Next we expand the arc length standing for the reduced 
surface energy. It is expressible as 



1 



Vl + C'(0 2 d£-1 



The most difficult part is the elastic energy. We start 
with the setup of the Airy stress function. Again, it 
is sufficient to go to order 2 in the ansatz known from 
Eq. (|l6|), which we have slightly modified in order to 
facilitate later calculations, 



9 + ~~~T Z ) cos(nkx)e nkz 
n*k z nk 



(A5) 
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APPENDIX A: EXAMPLE FOR THE 
CALCULATION OF THE ENERGY TERMS TO 
FOURTH ORDER IN THE AMPLITUDES 

In order to provide a taste of the algorithm, we present 
the method of calculating the necessary energy terms up 
to order four now. Basically this is a shorter (in terms of 
writing effort) but generalized version of EJ. Presuming 
a n r-j e™, using the reduced amplitudes A n and carrying 
out the transformation xk — ► £ we may write 



C(0 = * _1 I> n A.cos(nO 



(Al) 



and use the parameter e for bookkeeping. 013 and 04 will 
not play any role for the order 4, because any combina- 
tions of 03, 014, as well as a^ >n , a^n, ba <n , b^ n up to order 
4 in the energy terms disappear due to translational in- 



First, we write down the components of the additional 
stress tensor originating from the derivatives of x, m do- 
ing so assuming the orders of a n and b n to be e n . Addi- 
tionally we carry out the usual transformations kx — > x 
and kz — » z: 



n=l 



' cos (nx) (a n + 2b n + b n zn) e r ' 



Oxz = ^e"sin(nx) (a n + b n + b n zn) e nz , (A6) 

n=l 

2 

<j zz = — ^ e" cos (nx) (a n + b n zn) e nz . 

We now have a scheme for the definite calculation of as 
many Airy coefficients a^-, bij as necessary. We have to 
keep in mind that, according to our ansatz, the sum of 
subscripts of an Airy coefficient is identical with its order. 
For instance, 123,0 and ai,2 are both of order e 3 . 

Our scheme yields a set of equations at every order e" 
allowing us to calculate the coefficients whose subscript 
sum is n, respectively. The following steps have to be 
carried out: 



1. Write down the equations ( 14) using the stress com- 
ponents (A6) and the transformed derivative 



2. Expand those equations up to order z 3 . (A z 4 term 
would produce at least e 4 terms in the equations. 
These are not necessary for the calculation of the 
energy to order e 4 as can be seen from table |.) 
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Replace z by the interface (Al) 



Expand the Airy coefficients a n and b n to suffi- 
ciently high order according to Eq. (|l7|). This 



n + m < 3, omitting odd values of m. 
5. Keep every term in the equations up to order e 3 . 

According to every order e n we get two equations ( |A7[ ) , 
in which we have replaced sin(nkx) and cos(nfca;) by S n 
and C„ for brevity. Moreover we have used the theorem 
to be proven in Appendix [b| by setting all coefficients 
whose second subscript is odd to zero. 



Order 


Mode 


e 6 


T 


/ 


e 5 


[ai,4, &i,4] 6 


/ 


e 4 






e 3 


[ai,2,6i,2] 4 


[«3,0; ^3,o]g 


e 2 




[O2,0, &2,o] 4 


e 1 


[ai,o, &i,o] 2 








(S 2 ,C 2 ) (Ss.Os) — » 



TABLE I: Schematic view of the expansion method. It shows, 
which pairs of coefficients come out of which mode and order. 
Subscripts denote the energy order they are necessary for. We 
already took into account that coefficients with an odd second 
coefficient are zero. 



(Ai + bifi + ai,o) Si 
— <ii,o Ci 

(2Ai&i >0 + AiOi,o + h,o + a 2 . + 2A 2 ) S 2 
(Aibifi + Aiaifl + a 2 ,o) C 2 



0-1,2 + h,2 + g&l.oA + g°l,oA + 2°1. 0^-2 + &1,0^2 + &2,o4l + T^.O^-l J Si 

+ (equation for 63^) S3 

^a ltQ A 2 + ^bi,oA 2 - ^a 2fi A 1 - i& 2 ,o^i - 01,2 - ^1,0^1 - |&i,o^i^ Ci 

+ (equation for 03,0) C3 



= 



= 



(A7) 



Now, the coefficients of every function S n and C n have 
to be zero separately, which finally provides a system 
of six equations for the coefficients ai.o 5 o\ 2 , CL 2t o, 61,01 
&i j2 and 62,0 ( a 3,o and 63,07 which we would get from the 
coefficients of 63 and C3 are not necessary for order 4). 
The solution is 



0-1,2 



01,0 = . 

- A 2 

— A i ■ 

AiA 2 . 



a 2 ,o 
1 



61.0 = -Ai , 
6 2 ,o = A\- 2A 2 



(A8) 



-A 



-A-iA 



l-n-2 



Finally, the elastic energy is calculated by taking the 
reduced elastic energy density 1 /2((l + a xx ) 2 +a 2 z + 2cr 2 , z ), 
integrating first over z from —00 to C, expanding the re- 
sult to sufficient order in £, replacing £ by £(x), integrat- 
ing over x and keeping everything up to order e 4 . The 
result is the first line of Eq. ( |l9| ) . 



APPENDIX B: PROOF FOR THE VANISHING 
OF ODD EXPANSION COEFFICIENTS IN THE 
AIRY CONSTANTS 

The basis of this consideration are equations ( |l4] ) 
which we write down here again. We will call the l.h.s. 
of these equations ti and i 2 . 



fi := - (1 + a xx ) Cj{x) + a x 

t 2 := -UXZ C'(-T) + 0-22 . 



(Bl) 



The principle we have to apply is the generalized version 
of the one exhibi ted in a ppe ndix (|a|): Use the terms from 
Eqns. ( A2) and ( |A6| ) in (Bl), expand the Airy coefficients 
a n and b n as well as the exponential functions, replace 
z by the surface (Al), and then first sort by orders of 
e and then by modes. This will provide a system for all 
expansion coefficients of a n and b n . Note that in this cal- 
culation no square root is included. So the elastic field 
calculated in this way will be correct for all £, regardless 
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of the value of the derivative, for which the power series 
converges. Since the functions to be expanded are expo- 
nentials, for which the radius of convergence is infinite, it 
is expected that the expansion of the elastic fields has a 
much larger radius of convergence than that of the square 
root. 

Definition B.l The n-th order part of a term t will be 
denoted by v- n '. The selection of a mode m from such a 
term will be called I . 



Obviously, t\ and t% do not contain zeroth order parts. 



So, 



n-l 



tT ] =-£ 4™ ] C'(*) [,1 - ml + 41 - C'M W , (B2a) 



711 — 1 



4 n] = -£ 



(B2b) 



Consider any one of the stress terms. First take a look at 
the exponential function. In the course of our algorithm 
it will be expanded, i.e., every z 3 with j > will appear 
in the sum. Provided £(£) is given as in (Al), we note 
that C(a;)[ 2m ] corresponds to even modes and C,(x)^ 2m+1 ^ 

to odd modes. Consequently, when constructing C(£) 2 '^ 
using the well-known formula 



cos (nx) cos (mi) 



— [cos ((n — m)x) + cos ((n + m)x)] 



any even (odd) order j contains only even (odd) modes 
again. We can conclude by induction that this holds for 
any C(£) j an< ^ hence we formulate 
Lemma B.l Provided £(£) is given in the form of 
Eq. (Al), the expansion of the term exp(nC(£))^ ™^ 



contain only even (odd) modes, if j is even (odd). 

Without proof we additionally state 

Lemma B.2 exp(n((^))' : '' contains the cosine modes 
0-..J. 

For shorthand notations we now define 

exp(nC(0) [j1 
cos (n£) 
sin (n£) 



— K 

— C n , 
= S n . 



Next we consider an e-order of a generalized stress-like 
term, call it crfy . Let it contain some general Airy coef- 
ficient gj, so that 



°> = E 1 

i=i 

Of course, for the general form, we could have used sines 
instead of cosines as well. We expand the Airy coefficients 
and the exponential functions: 



9j = E "''o - > 

i=0 



/ , 3r 

i=0 



In] 

Now, o\lv has to be composed from parts of the n lower- 



most summands in this series: 

n n—j 



^ = Y. C iY,^n-i- j E j , i , n 1 ... x . (B3) 



It is useful to subdivide the terms (B3) into those with 
even and odd n: 

2n— 1 2n-l-j 

4 2 "~ 1] = E ^ E ffi,2n-i-*-A . (B4a) 

3=1 «=0 



2n 2n-j 

< n] = E ^ E g^-i-jEj.i . (B4b) 

3=1 t=0 

To be able to give brief descriptions we introduce 

Definition B.2 Airy coefficients whose second index is 
even (odd) will be called even (odd). 

By inspection we can derive now 



Lemma B.3 crjf" ^^z and cr|, 2 "'|2/-i contain odd Airy 
coefficients only, crjj 2 ™ ^laj— l and erj 2 ™' 1 2Z contain even 
Airy coefficients only. 



Lemma B.4 Lemma \B.H\ also holds for terms of the ex- 
tended form oiT ] C'(£) [2Tl-1-m] and <r [ $C0 2n ~ m] , re- 
spectively. 

Consequently, we can subdivide our set of equations 
for the Airy coefficients {t^ | m , \m} m to 

„ ._ r>[2n-l], .[2n]i J2n-1], ,[2n], , 

SX •— \Ii |2m,Ci |2m— lj^ |2m)t2 |2m— 1/ j 



•S2 



r [2n], ,[2n-l], [2n] , .[2n-l] , , 

ill 12m! 1 ! |2m— lj^a |2-roi r 2 |2m— 1/ 



The set si contains the equations which determine the 
odd Airy coefficients and S2 stands for the even ones. 

The only source of terms that are not products includ- 
ing Airy coefficients is the single £' (£)'"' hi Eq. ( B2a ). 
Now, C'(0^ i s simply equal to — nA„ sinn£, and conse- 
quently C0 n] \ m = -nA n S n , m . 

So, this term does not contribute to the ii-equations 
in s\ but only in s%. Hence si is a homogenous system 
for the determination of {a n ,2m-i, &n,2m-i}- It is solved 
by {a n> 2m-i = 0, & n ,2m— l = 0}, which had to be proven. 

Remark: It is possible but lengthy to show that this 
is the only solution. In principle one takes advantage 
of the fact that the determination of the odd coefficients 
can be carried out uniquely by a bottom-up scheme using 
the empty entries in table || 
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